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Abstract 

As galaxy redshift surveys probe deeper into the universe, they uncover ever 
more dramatic structures in the large-scale distribution of galaxies. In particular, 
the CfA2 and SSRS2 surveys to an apparent magnitude limit of 15.5 exhibit an 
impressive complex of sheets, filaments, and clusters. The formation of the large- 
scale structure in the universe results from the gravitational amplification of the 
primordial small perturbations of density. The primordial density perturbations are 
thought to be random fields originated as quantum fluctuations at the very early 
stage. Thus the understanding of the formation of the large-scale structure may 
reveal important information about the early universe and the laws of fundamental 
physics. One of the major obstacles to understanding the formation of the large- 
scale structure is the complexity of the evolution of the density inhomogeneities at 
the nonlinear stage when the observable structures form. One way of addressing this 
problem is to run three-dimensional numerical simulations. Here we review another 
approach based on the approximate analytic model of the nonlinear gravitational 
instability utilizing Burgers’ equation of the nonlinear diffusion. 


1 Introduction 

The term large-scale structure in the universe is referred to a distribution of galaxies on 
the scales roughly from I Mpc to fOO Mpc [17], where I Mpc = 10® pc ^ 3 ■ 10^“^ cm 
is a unit of length commonly used in cosmology. Galaxies can not probe much smaller 
scales because of discreteness, and on larger scales the galaxy distribution becomes almost 
homogeneous. The redshift surveys reveal spectacular abundance of structures often de¬ 
scribed as hlamentary, network, or bubble structure [17], [6], [4] (see Fig.I). The origin 
of the large-scale structure is one of the most important problems in modern cosmology. 
Many fundamental issues in physics, cosmology and astronomy ranging from speculations 
on the physical nature of dark matter, to the measurement of angular anisotropy of the 
microwave background radiation and determination of the epoch of galaxy formation join 
together here [21]. 
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The most popular and best developed class of theories of structure formation is based 
on the assumption that it started from primeval small amplitude density perturbations 
which grew by gravitational instability. Primeval perturbations are assumed to arise 
as vacuum fluctuations during the very early stage when the universe was expanding 
exponentially (inflationary universe) [13]. Afterwards, the density perturbations had a 
long history before they become galaxies, clusters of galaxies, superclusters and voids. The 
formation of galaxies is a very difficult problem itself. Many complex physical processes 
like star formation and supernova explosions are very important for understanding the 
galaxy formation. We shall discuss the mass distribution assuming that galaxies are fairly 
good (though not perfect) tracers of mass on large scales. 

As long as the density perturbations are small by amplitude their evolution is described 
by the linear theory of gravitational instability (see e.g. [28]; [18]; [22]. The linear theory 
is very well understood; in particular, it predicts the rate of growth of the perturbations 
at different stages of the evolution of the universe. Unfortunately, it is not applicable to 
the nonlinear stage when the amplitude of the density perturbations becomes large and 
the observable structures (sheets, hlaments, and clusters of galaxies) form. 

One of the most difficult obstacles arising in understanding models is the analysis 
of the evolution of perturbations at the nonlinear stage, when the typical amplitudes of 
density inhomogeneities become larger than mean density in the universe: hp/p > 1. 
The most straightforward way to address the problem of the nonlinear evolution is to 
do gravitational three-dimensional N-body simulations. Usually in simulations of this 
type the medium is assumed to consist of collisionless particles, in agreement with the 
popular hypothesis that the most of the mass in the universe is in the form of weakly 
interacting particles such as massive neutrinos or axions. The trajectory of each particle 
is calculated in the gravitational held generated by all the particles. Boundary conditions 
are commonly assumed to be periodic. 

Here we review an alternative approach to the problem of the large-scale structure 
in the universe. We present an approximate analytic technique to solve the nonlinear 
gravitational instability based on Burgers’ equation. Other analytic and semi-analytic 
methods are discussed in an excellent review by Sahni and Coles [19]. 

The initial condition is the result of the linear theory applied to the earlier stages. 
In the linear regime the density huctuations are assumed to be a Gaussian random held 
specihed by the spectrum and the amplitude. The current measurements of the angular 
huctuations in the temperature of the microwave background radiation by COBE (Cosmic 
Background Explorer) and other experiments put strong constraints on the initial huctua¬ 
tions. In particular, the amplitude of the temperature huctuations suggests that the scale 
where the density huctuations have recently reached nonlinearity is about 5 — 10 Mpc 
which is in a good agreement with the observations of the large-scale galaxy distribution 
(see e.g. [26]). 

According to the prediction of cosmic inhation, the geometry of the universe is assumed 
to be hat (see e.g.[13]. Eor simplicity we also assume that the cosmological term A = 0. 
Such a model is specihed by only two parameters: the Hubble constant Hq describing 
the current rate of expansion of the universe and the fractional density of baryons ftf,. 
The astronomical measurement of the Hubble constant is not very accurate: 50 < Hq < 
100 km Mpc~^] this uncertainty is usually expressed in terms of a parameter h\ 
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Hq = 100 h km Mpc~^ with 0.5 < h < 1. The density of baryons as well as the 
density of other components is convenient to measure in the units of the critical density: 

= pb!Per^ where the critical density per = SHqISttG ~ 2 • cm~^ is a parameter 

separating the closed cosmological models ftfot > 1 from the open ones ftfot < 1- An 
open universe has negative spatial curvature and expands forever, and a closed universe 
has positive curvature and eventually collapses. If ftfot = 1 the spatial geometry of the 
universe is flat and it expands forever. 

We shall assume that about 90% or more of mass in the universe is in a form of weakly 
interacting collisionless particles (dark matter) and the remaining few percent of the mass 
is in baryons; ftfot = ^dm + = 1 [21]. Although all luminous objects (e.g. stars) are 

made from baryons the baryon component dynamically is not very important on large 
scales. Therefore we shall study the evolution of density fluctuations in the collisionless 
medium neglecting the baryonic component. 


2 Basic Equations 


The evolution of density inhomogeneities can be described by a system of three partial 
differential equations comprising the continuity, Euler, and the Poisson equation (see e.g. 
[28], [18], [22]). It is convenient to exclude the homogeneous expansion of the universe 
from the analysis by introducing new variables. Instead of the coordinates r the comoving 
coordinates x = r jait) are commonly used, here a{t) is the scale factor describing the 
homogeneous expansion of the universe; in a flat universe a{t) oc Instead of the 

velocity u the so called peculiar velocity Up = u —a/a-r is used; if Up = 0 the velocity held 
on the scales smaller than the cosmological horizon {Rh ~ 3, 000 h~^ Mpe) is described 
by the Hubble law: u = ^ ■ r = H{t) ■ r). The velocities and gravitational held in the 
process we discuss are nonrelativistic [v c, p G) thus the evolution of density 
inhomogeneities can be described by the equations of classic hydrodynamics and gravity. 
In terms of the comoving coordinates and peculiar velocities the equations are as follows: 
the continuity equation 


^ , i' 

dt a 


the Euler equation 

and the Poisson equation 


9up 1 
—- + -( 
dt a 


a 

: -3-/9, 
a 

(1) 

-V</- Up, 

(2) 

a a 

- d), 

(3) 


where a = da/dt] p and p are respectively the density and mean density of mass; </> is the 
perturbation of the gravitational potential due to the inhomogeneities of density; G is the 
gravitational constant. 

The pressure is neglected since we study the medium interacting only gravitationally. 
Additional terms on the right hand side of the continuity and Euler equations are due to 
the homogeneous expansion of the universe and the factor 1/a is due to differentiation 
with respect to the comoving coordinates x: V = d/dxi = a ■ d/dvi. We assume that 
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the initial condition is small density and smooth velocity perturbations imposed on a 
homogeneous density distribution. 

As long as the amplitude of the density perturbations remains small their evolution 
can be analyzed in the linear approximation obtained by the linearization of the above 
equations. The exact solution of the linearized system has one growing mode which is 
the major object of our analysis. The velocity in the growing mode is a potential vector 
held. In the linear regime the spatial structure of the perturbations (in the comoving 
coordinates) remains unchanged and its amplitude is proportional to the growing solution 
bg of the differential equation 

d% db ^ 

a— + 2a— + 3a6 = 0. (4) 

at 

The scale factor a is assumed to be a known function of time and in a hat matter dominated 
universe bg[t) oc a[t) oc It is convenient to make yet another transformation of the 

variables: 


h(x,&p) 

= a p(x,t). 

(5) 

, , dx. 

1 ^ ^ 


v(X’6p) = — 

dbg 

= ^Vp(x,t), 

abg 

(6) 


= {3dabg)~^ (/)(x, t). 

(7) 


and also to reparametrize the time coordinate by the monotonic function of time bg[t) 
describing the growth of the perturbations. Finally, after explicit use of the function 
bg = a (which assumes ftfot = 1) and introducing the velocity potential $: v = — the 
equations take the form 

= 0 , ( 8 ) 

= ( 9 ) 

1 rj — fj 

a fj 

In the linear regime both the gravitational potential Lp and velocity potential $ remain 
constant and equal to each other, and the right hand side of eq.9 vanishes. In 1970 
Zel’dovich [27] suggested to extrapolate this condition well into the nonlinear regime; the 
corresponding solution is known in cosmology as the Zel’dovich approximation. 



^ + V.(,v) 
dv 

^ + ,v.V,v 

vV 


3 Zel’dovich Approximation 

The Zel’dovich approximation is convenient to formulate as a mapping from the La- 
grangian space L{q} into Eulerian space E{x} 

x(q, a) = q +a • Vo(q), ( 11 ) 

which obviously follows from eq.9 assuming that (/? sa $. The initial velocity potential 
$o(q) (t^ 08 '(q) = —d^o/dqi) is assumed to be a smooth random Gaussian held specihed 


4 



by the spectrum P<^{k). In cosmology the initial condition is usually characterized by 
the spectrum Ps{k) of the linear density fluctuations S = [p — p) j p = [r] — f])/f] which is 
obviously related to the spectrum of the potential 

P^{k) = k-^ Ps{k). (12) 

Utilizing the law of mass conservation one hnds the density as a function of the time 
coordinate a and the Lagrangian coordinate q 

fj 

i)(q, “)=[!_ „ (13) 

where Ai(q), A2(q) and A3(q) are the eigenvalues of the deformation tensor di^ = dqidq^ 

Combining eqs.ll and 13 one can hnd the density distribution in the Eulerian space. It 
follows from eq.ll and eq.l3 that the hrst objects form at the maxima of the largest 
eigenvalue and have very oblate shapes. After Zel’dovich they are known in cosmology as 
“pancakes”. Recent three-dimensional gravitational N-body simulations are in a perfect 
agreement with this conclusion [23]. The pancakes originate as the three-stream flow re¬ 
gions bounded by caustics, the surfaces of formally inhnite density. The shape and other 
characteristics of the pancakes are determined by catastrophe theory [1]. At the later 
stages the Zel’dovich solution predicts several different types of singularities which are 
classihed in [1]. 

We use this opportunity to remark that the Zel’dovich approximation (in two-dimensional 
space) is very similar to the equations describing the propagation of light in geometric 
optics [29]. 

The Zel’dovich approximation proved to be very good until orbit crossing when caustics 
form and the multi-stream flows occur (see e.g. [22] and references therein). (At this stage 
the original equations must be obviously modihed in order to incorporate the multi-stream 
flows.) However, the Zel’dovich approximation predicts the multi-stream flow regions to 
broaden very fast which contradicts to the results of the N-body simulations [5], [9]. 
Numerical studies of the orbits in the multi-stream flow regions show that the velocity 
component orthogonal to the pancakes randomizes very quickly [12]. As a result the 
pancakes observed in the N-body simulations remain quite thin. This result has become 
a physical basis for the adhesion approximation. 


4 Adhesion Model 

The general idea of the adhesion model is very simple. We wish to use the Zel’dovich 
solution everywhere except the regions of multi-stream flows. By adding a diffusion term 
into the Euler equation one can suppress the formation of the multi-stream flow regions. 
Assuming that the gravitational potential is approximately equal to the velocity potential 
(/? sa $ and adding a viscosity term in eq.9 one obtains the equation of the nonlinear 

diffusion [7], [8] 

—h (v • V)v = I/V^v. (14) 

da 
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Generally speaking the viscosity term need not to be in the form of eq.l4 but choosing this 
particular form one obtains Burgers’ equation that has an exact analytic solution [3]. For 
potential motion v = — eq.l4 can be solved by performing the Hopf-Cole substitution 
$(x,a) = —2i/log f/(x, a). As a result eq.l4 translates into the familiar linear diffusion 
equation 

f = ( 15 ) 

Solving eq.l5 for the velocity we obtain 


exp[-A(x,a;q)/2t/] 
/ d^q exp[—S'(x, a; q)/2i/] 


(16) 


where the “action” 

A(x, a; q) = -$o(q) + • (17) 

In cosmology the adhesion model has been used in two forms: one assumes a small but 
hnite value of the viscosity parameter u and the other assumes it is inhnitesimal: i/ —> 0 

[19]. 

For hnite u the trajectory of a particle can be determined by solving the integral 
equation [24], [25], [16], [15] 


x(q,a) = q + 


0 


da v[x(q, a^), a^]. 


and the resulting density can be determined from the continuity equation 


? 7 (x, a) = qjdet{ 



(18) 


(19) 


For an inhnitesimal value of the viscosity parameter i/ —> 0, the integrals in eq.l6 can 
be evaluated using the method of steepest descents [7], [8], [10], [11]. In this case 


v(x, a) 


X — q(x, a) 
a 


( 20 ) 


where q(x, a) is the coordinate of the absolute minimum of the action 5'(x, a; q) at given 
X and a. The points q that minimize the action obviously satisfy the Zel’dovich equation 
11 . 

This solution (eq.20) has an interesting geometrical interpretation (see Fig.2). One 
can hud the Fulerian coordinate x of the particle with initial Lagrangian coordinate q at 
a chosen time simply by projecting the apex of the paraboloid 

P(x.o;q) = (21) 

a 


assuming that there is a constant C satisfying simultaneously two conditions: 1) the 
paraboloid P is tangent to the initial velocity potential $o at q and 2) it does not cross 
$0 at any point. At early times a is small and the curvature of the paraboloid is greater 
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than that of $o and the above two conditions can be easily fulhlled. The Zel’dovich 
approximation is universally valid at this stage. 

As time passes the curvature of the paraboloid decreases. As a result the points q 
appear that do not satisfy the above two conditions. Such a point has been stuck into a 
surface due to orbit crossing. Again, such surfaces can be found by projecting the apex 
of the paraboloid but this time the paraboloid is tangent to the initial velocity potential 
in two points simultaneously. These surfaces correspond to the pancakes; their thickness 
is proportional to the value of the viscosity parameter and therefore is inhnitesimal as 
1 / —> 0 . 

Later two more types of structures form: the hlaments and knots. The hlaments 
correspond to the case when the paraboloid touches the initial potential in three points 
simultaneously and the knots when it touches the potential in four points. The set of 
surfaces, hlaments and knots make a cellular structure: large regions of low density are 
separated by the surfaces, the hlaments are at the intersections of the surfaces, and the 
knots are at the intersections of the hlaments. This geometrical construction can be 
viewed as the skeleton of the real structure. In one-dimensional case illustrated by Fig.2 
obviously there are no surfaces nor hlaments. 


5 Accuracy of the Adhesion Model 

Apart from the regions of high density where the viscosity term (eq.l4) plays a signihcant 
role the adhesion model is exact in one-dimensional case. It means that the velocity held 
outside the high density regions is predicted exactly. Also the motion of the clumps is 
described very accurately. This is because the Zel’dovich solution is exact outside the 
multi-stream how regions in one-dimensional case [22]. 

In more interesting two- and especially three-dimensional case the adhesion model is 
only an approximation. Therefore the question arises about the accuracy of the model. 
The both variants of the adhesion model have been thoroughly tested against the grav¬ 
itational N-body simulations in two and three dimensions. Both the N-body simulation 
and the adhesion model used the identical initial conditions and were compared at several 
stages of the evolution. 

The geometrical version of the adhesion model was tested against the two-dimensional 
N-body simulations with the initial power law spectra Ps{k) oc A;” with spectral indices 
n = 2,0 and —2 and various cutoffs [11] 

The N-body simulations used the particle-mesh code with 512^ particles on equal 
mesh and periodic boundary conditions (for the details see [2]). The code constructing the 
skeleton of the structure is described in [10]. It has been found that the skeleton reproduces 
the density distribution extremely well (see Fig.3) for all choices for the parameters of 
the initial spectra until the stage when the scale of the nonlinearity A;”/ dehned by the 
equation 

rkni 

/ Ps{k) d^k = 1 (22) 

Jo 

reaches the characteristic scale of the initial velocity potential $o: k~i^ < Rq>g, here D is 
the dimensions of the space. The scale of the potential is dehned from the expansion of 
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( 23 ) 


the correlation function +•••)• 

Rt, = {2D)'I‘ 

c^l 

where (Tq and (Ti are the dispersions of the potential and its gradient, respectively. In 
the N-body simulations we deal with hnite ranges and therefore always exists. At 

the later stages k~l > the adhesion model remain qualitatively correct though its 
accuracy somewhat deteriorates. 

The version of the adhesion model utilizing a hnite viscosity parameter u has been 
quantitatively compared to fully nonlinear, numerical three-dimensional gravitational N- 
body simulations [15]. The initial perturbations were again random Gaussian holds with 
power-law spectra with the spectral indices n = —2, —1, 0,-|-1. The particle-mesh N- 
body and adhesion simulations both used 128^ particles on a 128^ mesh and periodic 
boundary conditions. In these simulations of the adhesion model, the smallest value of the 
viscosity parameter u that did not produce numerical overhows has been used. For further 
discussion of the N-body and adhesion simulations see [14] and [25] respectively. The both 
codes calculated the particle positions therefore the corresponding density distributions 
could be easily generated (Fig.4). The primary tool of the quantitative comparison was 
the cross-correlation coefficient for the density helds obtained from the gravitational N- 
body simulation and the adhesion model with the identical initial conditions. Also the 
density distribution functions and the power spectra were compared. The adhesion model 
produces an excessively hlamentary distribution due to smoothing effects in high density 
regions. As a result the density distribution function in the adhesion model is lower than 
that of the N-body simulation at high densities pjp > 10 — 20 depending on the initial 
spectrum and the power spectrum of the nonlinear distribution falls off steeper on small 
scales k > kni- 


6 Summary 

The adhesion model based on three-dimensional Burgers’ equation of the nonlinear diffu¬ 
sion has been found to work very well in explaining the large-scale features of the structure 
of the universe. Comparison with gravitational two- and three-dimensional N-body simu¬ 
lations has shown remarkable agreement till very late nonlinear times. The main drawback 
of the adhesion approximation is probably the fact that it cannot describe accurately the 
density distribution within pancakes and hlaments. 

The adhesion model provides a natural qualitative explanation of the origin of the 
large-scale coherent structures such as superpancakes and superhlaments, as a result of 
coherent motion of clumps due to large-scale inhomogeneities in the initial gravitational 
potential (compare Fig.l and 5). The formation and evolution of large-scale structure is 
described by the adhesion model as a two stage process [20]. During the hrst stage matter 
falls into pancakes and then moves along them towards hlaments and then along hlaments 
to collect hnally in knots. At the end of the hrst stage the formation of the skeleton of 
the large scale structure is complete and virtually all of the matter in the universe is 
located in one of three structural units: pancakes, hlaments or knots. The second stage 
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sees the deformation of the large-scale structure skeleton due to the dynamical motion of 
pancakes, hlaments and especially knots. At this stage knots merge into larger knots and 
small voids disappear giving space to growth of the larger ones. Eventually almost all 
the mass concentrates in knots. Depending on the initial spectrum the knots may move 
coherently in such a manner that they concentrate to superpancakes and superhlaments. 
The superpancakes and superhlaments can be identihed by applying the adhesion model 
to smoothed initial potential [11]. 
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Figure 1: This is a low resolution version of Fig.2 from [4], It displays galaxies from 
two redshift catalogs to an apparent magnitude limit of 15.5 and distances less than 
120 h~^ Mpc\ CfA2 (north; upper portion) and SSRS2 (south; lower portion). For CfA2 
the box shows the declination limits 8°.5 < 8 < 44°.5 and the right ascension limits 
8^ < a < 17^. In the south the box limits are —40° < 8 < —2°.5 and the right ascension 
limits are 20^.8 < a < 4^. There are 9325 galaxies in the original image. 
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Figure 2: Geometrical prescription of descending a paraboloid onto the initial velocity 
potential in order to find the Eulerian positions of particles and knots in one-dimensional 
case. The particle having Lagrangian coordinate go in the uppermost panel has the 
Eulerian coordinate of the paraboloid apex. In the middle panel corresponding to a later 
stage it is stuck into the knot mi which velocity is shown by the arrow. The knot m 2 
is determined by the second paraboloid in the middle panel. The lower panel shows the 
knot M formed as a result of merging mi and m 2 . Adapted from [20]. 
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Figure 3: Composite picture of the results of the two-dimensional N-body simulation and 
the adhesion model for the flat initial spectrum Ps{k) oc A;°. The left hand side panel 
shows the model with the cutoff of the initial spectrum at kc = 32 kj^ here kj is the 
fundamental frequency corresponding to the box size: L^ox = 2tt jkj] the right hand side 
panel shows the model without a cutoff (except at the Nyquist frequency k^y = 256). Not 
all the particles can be shown, but this is a fair sample). Solid lines and circles represent 
the skeleton of the structure constructed by the paraboloid technique. The area of circles 
is proportional of the mass of knots. Adapted from [11]. 
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Figure 4: A gray scale plot of thin (iv6oa;/128) slices through the simulation cubes for 
n = 0 (left hand side panels) and n = —1 (right hand side panels) initial spectra, at 
the stages when kni = Skj. The gravitational three-dimensional N-body simulations are 
shown in the bottom panels and the adhesion model simulations are shown in the top 
panels. Adapted from [15]. 
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Figure 5: Simulated galaxy distributions are drawn from the adhesion model simulation 
based on the biased Cold Dark Matter cosmological model. These redshift-angle projec¬ 
tions show all “galaxies” with apparent magnitude less than 16.5 and distance less than 
180 h~^Mpc. Adapted from [24]. 
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